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Abstract. Plerionic supernova remnants exhibit radio emission with remarkably flat spectral indices ranging from 
q = 0.0 to a = —0.3. The origin of very hard particle energy distributions still awaits an explanation, since shock 
waves generate particle distributions with synchrotron spectra characterized by a < —0.5. Acceleration of high 
energy leptons in magnetohydrodynamic turbulence instead may be responsible for the observed hard spectra. 
This process is studied by means of relativistic test particle calculations using electromagnetic fields produced 
by three-dimensional simulations of resistive magnetohydrodynamical turbulence. The particles receive power-law 
energy spectra Nfa) oc with s ranging from 1.2 to 1.6, i.e. particle spectra that are required to explain the 
radio emission of plerions. 
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1. Introduction 

Fillcd-center supernova remnants (SNR), or plerions ex- 
hibit flat radio spectra, with power-law indices < a < 
0.3 (fo r S v o c v~ a ) (Weiler & Panagia EZHI Weiler & 
Shaver 1978). Such spectra require energy distributions 
of the electrons Nfy) oc 7 _s with s = 1 + 2a be in the 
range 1 < s < 1.6. The flatness of their radio spectra dis- 
tinguishes plerions from a typical shell-type SNR with a 
mean a ~ 0.5 implying AT (7) oc j~ 2 (e.g. Green 1991). 
The origin of relativistic electrons responsible for these 
unusual radio spectra of plerions is not yet understood 
(e.g. GreenHHa Woltjer et al. lTW71 Arons lP39g|l . 

Although the many details of the nature and structure 
of plerions are still unclear (see Arons 119981 Salvati et 
al. 1998. Bandiera 2002 for reviews) there is a consensus 
on the following scenario: Plerions are expanding bubbles, 
filled with magnetic fields and relativistic leptons. Both 
components are continuously supplied to the nebula by 
some central source, i.e. by a rotating neutron star in form 
of a strongly magnetized wind whose energy is at least 
partially dissipated in termination shock waves (Kennel 
& Coroniti lT5H2alll984bl Galant & Arons DHH) . 

Arons (1998) discussed the difficulties of particle accel- 
eration in plerions in great detail. With respect to proto- 
typic Crab nebula he pointed out that the emission at high 
frequencies (from the optical towards X- and 7-rays) di- 
agnoses the coupling physics today since the synchrotron 
loss times of the high energy particles are significantly 
smaller than the lifetime of the nebula. The radio emis- 
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sion instead measures the integral of the pulsars input 
over the entire history of the nebula - most of the stored 
relativistic energy is in magnetic fields and radio emit- 
ting particles (about 1Q 50 ergs). Arons noted that aver- 
aged over the whole nebula, the radio emitting spectrum 
has the form AT (7) oc 7- 15 , 10 2 5 < 7 < 10 4 based on the 
detailed spectral index maps by Bietenholz & Kronberg 
( 1992) which show the particle distribution to be remark- 
ably homogeneous. The origin for this energy distribution 
is not clear since the wind termination shock wave models 
do not yield power laws flatter than AT (7) oc j~ 2 at en- 
ergies small compared to 7 ~ 10 6 . Arons speculated that 
some additional acceleration physics by magnetohydrody- 
namic (MHD) waves may have to be included. 

Bandiera et al. lj2tl0 2) mapped the Crab Nebula at 230 
GHz and compared it 1.4 GHz map. The spectral index 
in the inner region at 230 GHz is flatter (by ~ 0.05) than 
in the rest of the nebula. Furthermore they found a steep- 
ening of the spectrum at the places of radio emitting fila- 
ments, concluding that the magnetic field strengths in the 
filaments is higher than in the surrounding nebula. Some 
evidence for in-situ acceleration in the Crab nebula has 
been collected by the radio observations of Bietenholtz, 
Frail & Hester (2001). 

It is the aim of our contribution to show that resistive 
magnetohydrodynamical fluctuations excited to a highly 
turbulent level accelerate test particles to energy distribu- 
tions with power law indices 1 < s < 1.6 being in good 
agreement with the radio observations of plerions. We also 
show that an increase in particle density as can expected 
for the filaments lead to a steepening of the particle spec- 
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trum. Given the fact, that central sources energize the neb- 
ula, thereby exciting strong MHD fluctuations we think 
that our simulations may shed some light on the probably 
necessary additional acceleration physics as speculated by 
Arons l(TM5jl . 

In the next section we describe the MHD simulations. 
Section |31 contains the results of relativistic particle sim- 
ulations including radiative synchrotron losses and finally 
we discuss our findings in Section 4 

2. Reconnective MHD turbulence 

Pulsars permanently power their environment by radia- 
tion and, in particular, pulsar plasma winds. In fact, in 
the case of the Crab nebula the permanent energy input 
by plasma flow is about 10 38 erg/s. This plasma flow en- 



counters the termination shock at radius R s 
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where it gets decelerated and leptons are accelerated to 
ultra-relativistic energies. Most of their synchrotron radi- 
ation is emitted in the optical to 7-ray band. The presence 
of this termination shock is accompanied by strong exci- 
tation of Alfven waves which are expected to excite strong 
MHD turbulence in the pulsar's nebula (Kronberg et al. 
I1993fl . In this picture the size of the turbulence cells in 
the nebula is limited by i? s . The projected magnetic field 
within the Crab nebula was found to have a coherence 
length of about 10 16 cm (Bietenholz & Kronberg 1992). 
Kronberg et al (1993) suggested that the fact that this 
length scale is close to the inner shock radius R s , may be 
more than coincidental. Laboratory studies of turbulent 
flows have revealed the existence of persistence structures, 
which are advected in the flow and exhibit sizes compa- 
rable to the source size (Hussein 119831 I1986[) . In other 
words, the source of turbulence within a supernova rem- 
nant is given by a shock wave, like in the Crab nebula, the 
spatial scale of the shock is the maximum length scale for 
turbulence within the nebula. To be precise, in the Crab 
nebula the characteristic size of the turbulence should be 
smaller than 10 17 cm. 

The associated MHD turbulence can be numerically 
generated by means of the so-called Orszag-Tang turbu- 
lence (Orszag & Tang ll97^)l . The Orszag-Tang initial con- 
dition is a generic way to excite turbulence in a magne- 
tized plasma and is given by the non-linear interaction of 
Alfven waves which are MHD eigenmodes of a magnetized 
fluid e.g., 



B x = -B sm[2ir(y - y min )/(y max - y min )] 

sin [2n(z - z min )/ (z max - z min )] 
B y = S sin[47r(a;-x m i n )/(x max -a; m i n )] 

sin [2n(z - z min )/ (z max - z min )] 
B z = B sin[47r(x-a; m i n )/(a; max -x min )] 

sin [2n(y - y m i n )/(y max yminyj 
v x = -v Q sm[2ir(y - y mhl )/(y max - y mhl )} 

sin [2n(z - z min )/ (z max - z min )} 
v y = v sin [2-7r(x - x min )/(x max - x min )] 



(1) 
(2) 
(3) 
(4) 
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The amplitudes are chosen as Bo = 5.5 x 10 3 G and 
v = 0.05, 0.016, 0.005, 0.0016, 0.0005, 0.00016 c for the 
six different test particle runs discussed in section [21 

The nonlinear interaction of the Alfvenic perturbations 
result in almost homogeneous turbulence. We use peri- 
odic boundary conditions in all directions and a numer- 
ical box given by 

•^mini l/min: ^min — 0-5 X 10 Cm and 

£maxi 2/max, z max = 0.5 x 10 16 cm with a resolution of 101 3 
grid points. 

We model this kind of turbulence by means of a 
well approved resistive compressible 3D MHD code (Otto 
1990}. It integrates the balance equations that govern the 
macroscopic low-frequency dynamics which read 

dt 
dpv 
~~dt 



V • (pv) = 
+ V • (pvv) = 



-Vp- 



i-(VxB)xB 



dp _ 
_ = _ v . V p- 



7pV ■ v + (7 - 1) T) (V X B)" 



1 dB 

c~dt 



1 c 
= V x - (v x B) - —V x (77V x B) 

C 47T 



(7) 
(8) 
(9) 
(10) 



where p, v, p, and B denote the mass density, bulk ve- 
locity, thermal pressure, and the magnetic field. By r\ the 
resistivity is denoted. 

Reconnection is allowed due to a microturbulent resis- 
tivity which is modeled in a current dependent way. The 
motivation for this kind of resistivity is the following: In 
the highly turbulent almost ideal plasma of the nebula 
collisional resistivity is negligible. Rather a wide variety 
of micro-instabilities driven by electric currents is respon- 
sible for localized dissipation. For strong turbulence the 
amplitude of the resistivity is given by r\ sa m e w e /ne 2 , 
where uj e — Aime 2 / m e is the electron plasma frequency 
(Vasyliunas 1975). This kind of resistivity is assumed for 
our simulations. 

As a result of the resistive MHD calculations we get 
a turbulent field configuration which we use for the test 
particle simulations. Fig. ^ an d 121 show the magnetic en- 
ergy density and the magnetic field lines. These figures are 
to illustrate the turbulent structure of the magnetic field. 
Fig.|3shows the reconnection regions, i.e. places where the 
parallel component of the electric field is appreciable high. 
These are the acceleration sites where we expect the par- 
ticles to gain energy. They are also distributed in a very 
stochastic way clearly indicating the turbulent nature of 
our fields. 

3. Relativistic particle simulations 

The MHD calculations described in the previous section 
are used as the electromagnetic environment for the stud- 
ies of electron acceleration by means of test particle simu- 
lations. We want to know how electric particles behave in 
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Fig. 1. 3D Isosurface plot of the magnetic energy density. 
The shaded volume surface corresponds to a magnetic en- 
ergy density of u ma gn ~ 2.5 x 10 _11 erg/cm 3 . The spatial 
extension is given in units of 10 16 cm. 




Fig. 3. Reconnection regions, i.e. regions with appreciable 
value of Eu . The shaded volume surface corresponds to a 
value of E\\ w 1.3x 1CT 6 , 4x 10~ 7 , 1.3x 1CT 7 , 4x 10~ 8 , 1.3x 
10 -8 , 4x 10 -9 statvolt/cm for the six different test particle 
runs. 



Fig. 2. Magnetic field lines. 

the complex three-dimensional turbulent electro-magnetic 
field configuration, in particular how their energy distri- 
bution develops. 

The electric field is derived from the MHD quantities 
B, v and rj by means of the normalized Ohm's law 

1 c 

E = — vxB + — nV x B, (11) 

C 47T 

and gives for the mean value of of the electric field 
E w 0.02 B. Equation Ijlll) states that a parallel com- 
ponent of the electric field En only occurs for r\ ^ 0. 
Compared to the mean perpendicular component E± the 



mean parallel component in our simulations is found to be 
E\\ re 0.04 E±. 

The results of the MHD calculation are scaled to model 
different physical conditions, i.e. different mass densities 
inside the nebula, and are then used for the test par- 
ticle simulations. Here we present the results of parti- 
cle simulations for six different mean mass densities p = 
10 n m p /cm 3 , n = —2,-1,0,1,2,3 (corresponding to the 
initial vq). The mean magnetic field strength for all runs 
was chosen to be B = 3 x 10" 4 G. 

Since the data for the electric and magnetic field is only 
available on a discrete three-dimensional grid the test par- 
ticle code uses linear interpolation to determine the field 
values for any location. We note that the term "test par- 
ticle" means that the electromagnetic fields produced by 
the particles are not changing the global fields, though the 
effect of synchrotron radiation on the motion of the par- 
ticles is taken into account, i.e. the relativistic equations 
of motion have the form 

= q (E + ^ P xB) + F Rad , * = JL, (12) 
at \ jmc ) at jm 

where m and q denote the mass and charge of the particles, 
which are electrons in our case and 

t = \A+(— ) 2 ( i3 ) 

V Vtoc/ 

is the Lorentz factor. 

The calculations are relativistic, including the energy 
losses via synchrotron radiation and inverse Compton 
scattering by (Landau & Lifshitz fl951(l 

FRad ~ 2 3^{ EXB + 1 BX{BXV) + 1 E{E - V) 
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The equations are numerically integrated by a Runge- 
Kutta algorithm of fourth order with an adaptive stepsize 
control. As a result we get the momentum and location at 
certain times of each particle in the given ensemble. The 
codes has been used before to study high-energy particle 
acceleration (Nodes et al. 2003, Schopper et al.[1999). 

The direction of the initial momentum vector is equally 
distributed in a cone of opening angle 80° at the bottom 
of the box. This initial condition corresponds to the idea 
that the electrons enter the the turbulent region moving 
in a well defined direction. The notion "final state" relates 
to a set of particles in which each individual particle has 
left the computational box. 

The simulations start off with electrons injected into 
the computational box at z = z m in- We use a powerlaw 
^kin ran ging from 10 2 to 10 4 m e c 2 as the initial energy 
spectrum for all calculations. The choice of this injection 
spectrum was motivated by our synchrotron model for the 
infrared to X-ray emission of the Crab pulsar (Crusius- 
Watzel, Kunzl & Lesch^UMJ. There we could show that 
the emission of the neutron star could be explained with 
a single energy distribution N(E) oc E^ in . We could also 
show that such an energy distribution is naturalley pro- 
duced in an efficient pair cascade. Given that plerions in 
general are powered by a pulsar like in the case of the 
Crab nebula, an energy distribution oc E^ 2 n seems a quite 
reasonable choice as a probable injection spectrum in tur- 
bulent plasmas of a plerionic supernova remnant. 

Another series of simulations we performed showed 
clearly that the resulting energy spectrum is independent 
of the index of the initial energy distribution. Also the 
initial energy range does not play a crucial role for the fi- 
nal distribution. After the initial injection of the electrons 
they get accelerated very rapidly and their energy dis- 
tribution reaches a stationary state. This happens within 
only a few days and then the distribution remains for sev- 
eral years until the final state is reached. For different 
initial distributions only this first acceleration phase gets 
shorter or longer but the final distribution is given by the 
energy offered by the MHD fields. Thus, even for a non- 
relativistic initial energy distribution we expect the final 
spectral index to be the same as the one found in the 
relativistic runs. 

Fig- El shows the initial and the final energy spectra of 
the test particle runs for different values of the mean mass 
density ~p. The calculations show that for all values of p 
the spectrum gets flatter than the initial spectrum which 
means that particles are accelerated. Compared to the ini- 
tial spectrum the low energetic electrons get redistributed 
to higher energies, i.e. the low energy peak of the initial 
spectrum is reduced. 

The cut-off energy in the spectrum is correlated with 
the maximum energy that the electrons can gain from 
the mean parallel component of the electric field E'^ x = 
elE\\ , where I is the length of the computational box. 



The final energy spectra can be approximated by a 
broken power law consisting of two parts: a low energy 
part which has a relatively flat gradient and a steeper high 
energy part. Basically the latter part of the spectrum is 
responsible for the synchrotron radiation therefore we will 
only consider this part. 

For the first run, the one with the lowest density ~p 
we get a spectral index of —1.19 for the high energy part 
above energies of 7 w 4.0 10 3 . The spectral index s de- 
creases with increasing p (approximately by a logarithmic 
relation), i.e. the spectrum gets steeper for denser regions. 
This result corresponds well with the observed steepening 
of the radio spectrum at places of magnetic filaments. 

The resulting synchrotron spectrum can be calculated 
as the sum of total power per frequency emitted by each 
individual electron. The latter is given by Rybicki & 
LightmanEHT! 

J (-) = ^-Bsine- r ff 6/3 (0de, (15) 
V v c ) mc 1 v c J v/Vc 

where v c — -^^B^ 2 sm8 is the critical frequency, 8 is 
the pitch angle and K$ is the modified Bessel function of 
order 5/3. Fig.[5]shows the synchrotron spectra of the high 
energy particles calculated for the "real" pitch angle dis- 
tribution taken from our simulation and a isotropic pitch 
angle distribution. The isotropic distribution leads to a 
synchrotron spectrum with spectral index a — (s — l)/2 
where s is the index of the energy spectrum. This means 
that the steepening of the energy spectrum with increas- 
ing density results in a steepening of the synchrotron spec- 
trum. 

The actual final pitch angle distributions in our simu- 
lations are clearly anisotropic. In the low density runs the 
pitch angle distributions offer a peak at sin (9 ss 0.01. Since 
J(v) ~ sin# the total radiation power emitted is lower as 
compared to the isotropic case (fig. EJ- The pitch angle 
is also correlated with the gained energy, i.e. the highest 
energetic particles have the lowest pitch angles. This cor- 
relation results from the acceleration along the magnetic 
field lines for which the parallel component of the elec- 
tric field E\\ = (77V x B)« is responsible. E\\ can only be 
found at places with significant value of 77. Contrary to 
second-order acceleration processes the acceleration due 
to magnetic reconnection is able to explain the observed 
spectra by particle acceleration in turbulence. To prove 
the efficiency of acceleration in presence of parallel eletric 
fields we performed a test particle run without any En . 

Fig. H3 shows the results of a run that is similar to 
the first run in Fig. 0| and Fig. [5] with the modification 
En =0. The resulting spectra are considerably steeper, 
since the acceleration is less efficient as in the case with 
reconnection. Therefore the presence of parallel electric 
fields is crucial to explain the observed flat spectra. This 
result is also corroborated by particle simulations starting 
from an analytic Kraichnan type ideal MHD turbulence 
(these results will be published elsewhere). 



Nodes et al.: Particle Acceleration In Plerionic Supernova Remnants 



energy spectrum 



energy spectrum 



s = -1.19 




10 10 10 10 10 

r-i 




10 10 10 10 10 

r-1 



(a) p = 0.01 m p /cm 3 
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(c) p = 1 mp/cm 3 



(d) p = 10 mp/cm 3 



energy spectrum 



energy spectrum 





(e) p = 100 mp/cm 3 



(f) p = 1000 m p /cm 3 



Fig. 4. The initial (dotted line) and final (solid line) energy spectrum of the test particle simulations for different 
values of the mass density (p = 10 n m p /cm 3 , n = —2,-1,0,1,2,3). A least squares fit of the spectrum in a certain 
energy range is plotted above the spectrum and is labeled with the spectral index s. 



Nodes et al.: Particle Acceleration In Plerionic Supernova Remnants 



synchrotron spectrum 



10 10 

V [Hz] 



a = -0.10 




synchrotron spectrum 




10 10 
V [Hz] 



(a) p = 0.01 nip/ cm 3 



(b) p = 0.1 ?7ip/cm 3 



synchrotron spectrum 



synchrotron spectrum 



a - -0.15 




v [Hz] 




(c) p = 1 mp/cm 3 



(d) p = 10 m p /c 
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(e) p = 100 wip/cm 3 



(f ) p = 1000 mp/cm 3 



Fig. 5. The initial (dotted line) and the final (solid line) synchrotron spectrum of all runs calculated as the sum of 
the individual synchrotron spectra of the corresponding electrons (equation l|15fl ). The dashed line depicts the final 
synchrotron spectrum of a electron population with isotropic pitch angle distribution. The power per unit frequency 
J(v) is normalized to the total rest mass energy Nmc 2 of the set of electrons, where N is the number of all electrons. 
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energy spectrum synchrotron spectrum 




Fig. 6. Particle energy and synchrotron spectra of a single run with En = 0. The value for the mean density in this 
run is the same as in the first run in Fig. 0] and Fig. El (p = 10~ 2 m p /cm 3 ). Due to the lower mean electric field energy 
the achieved energy levels are lower. 



4. Conclusions 

We studied the acceleration of relativistic electron popu- 
lations in three-dimensional reconnective turbulence with 
applications to the Crab nebula. By means of test parti- 
cle simulations performed within non- linearly evolved elec- 
tromagnetic fields modeled by resistive MHD simulations 
and including radiative losses via synchrotron emission we 
could show that particles are accelerated and form a en- 
ergy distribution with powerlaw index a ranging from 1.2 
to 1.6. We can give limits for the maximum size of the tur- 
bulent regions in the nebulae. If the volume is too large 
and the particle escape time is too long their spectra get 
significantly harder, i.e. the spectral index is flatter than 
1 (see Fig. [7J| . For the case of the Crab nebula we can rule 
out turbulent cells larger than 10 16 cm. This limit is in ac- 
cordance with experimental results from turbulent fluids, 
in which the turbulence size is limited by the spatial scale 
of the turbulence source. Since the Crab nebula is pow- 
ered by a pulsar wind which is shocked at about 10 17 cm 
the spatial scale of the turbulence should be smaller. 

We like to suppose the following scenario for the ra- 
dio emission of the Crab nebula: The nebula is perma- 
nently powered by the plasma wind originating from the 
pulsar. Consequently strong MHD Turbulence is exited 
in the nebula. In the turbulent regions the leptons expe- 
rience in situ acceleration in the reconnective turbulent 
environment. Our simulations prove that the resulting en- 
ergy spectra and synchrotron spectra are considerably flat 
and therefore can explain the observed hard radio spectra 
(Wciler & Shaver H978[) . Whereas we concentrate on the 
Crab nebula, our findings should also be applicable to flat 
radio spectra of plerions in general. We note, that the effi- 
cency of the acceleration within th test particle approach 
depends on the ratio of Alfven speed to velocity of light 
which determines the effective electric field strength re- 
sponsible for particle energization. If the Alfven speed de- 
creases since the magnetic field decreases the electric field 
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Fig. 7. The initial (dotted line) and final (solid line) en- 
ergy spectrum of a test particle run with larger box size 
of I = 10 17 cm. The density for this run was chosen to be 
p = I0~ 2 m p /cm 3 . 



will also be weaker, reducing the acceleration to higher 
energies. If a supernova remnant consists of region with 
significantly varying magnetic fields and densities the non- 
thermal particles may be a superposition of particles ex- 
periencing different efficencies and electromagnetic fields. 

We note that Atoyan l|1999f) considered the origin of ra- 
dio emitting electrons in the Crab nebula in terms of inten- 
sive radiative and adiabatic cooling in the past. Whereas 
he argued that effective in situ acceleration in the nebula 
is not necessary to explain the radio spectra, our simu- 
lations prove the capability of magnetic reconnection to 
energize an electron population of a magnetized plasma 
under the influence of external electromagnetic fields onto 
which the accelerated leptons have no back reaction. 
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